/*=============================================================

PROGRAMMING TECHNIQUES FOR SUPERCOMPUTERS

Assignment Number 1

Task 2

=============================================================*/

#include <stdio.h>
#include "timing.hh"
#include <iostream>
#include <math.h>
#include <fstream>

int main(int argc, char *argv[]) 
{
	double wcTimeStart,wcTimeEnd,cpuTimeStart,cpuTimeEnd,wcTime, cpTime;	
	
	// Plot A for Vector Triad
	std::ofstream fileouta("plota.dat");
	fileouta << "# Plot for Task (a)" << std::endl;
	fileouta << "# N\t MFLOP/sec" << std::endl;

	// Plot B for constant Multiplication
	
	std::ofstream fileoutb("plotb.dat");
	fileoutb << "# Plot for Task (b)" << std::endl;
	fileoutb << "# N\t MFLOP/sec" << std::endl;
	
	// Plot C for Matrix Vector Multiplication
	std::ofstream fileoutc("plotc.dat");
	fileoutb << "# Plot for Task (c)" << std::endl;
	fileoutb << "# N\t MFLOP/sec" << std::endl;


	const double s = 2.0;
	std::cout << "Running Task 2.1 and 2.2" << std::endl;
	for (int k =3; k<21; ++k)
	{
		int N = int(pow(2.5,k));
		double *a = new double[N];
		double *b = new double[N];
		double *c = new double[N];
		double *d = new double[N];
		
		for (int i = 0; i<N; ++i)
		{	
			b[i] = i*2.0;
			c[i] = i*3.0;
			d[i] = i*1.0;
		}
	
		// Task 2.1: Vector Triad a[i] = b[i] + c[i]*d[i]	
		int repeat =1;
		double runtime =0.;
		for (; runtime <0.1; repeat*=2)
		{	
			timing(wcTimeStart,cpuTimeStart);
			for (int r=0; r<repeat; ++r)
			{
				for (int i = 0; i<N; ++i)
				{
					a[i] = b[i] + c[i]*d[i];
				}
				if (a[N-1] <0.) std::cout << "Dummy Loop Encountered"<< std::endl;
			}
				timing(wcTimeEnd,cpuTimeEnd);	
				runtime=wcTimeEnd-wcTimeStart;				
		}
		repeat/=2;
		long noFlops = N*2*repeat;
		double mFlops = noFlops/runtime/1000000;
		fileouta << N <<"\t" << mFlops << std::endl;
		
		// Task 2.2: Vector update: a[i] = s * a[i] // s scalar
		repeat =1;
		runtime =0.;
		for (; runtime <0.1; repeat*=2)
		{	
			timing(wcTimeStart,cpuTimeStart);
			for (int r=0; r<repeat; ++r)
			{
				for (int i = 0; i<N; ++i)
				{
					a[i] = a[i]*s;
				}
				if (a[N-1] <0.) std::cout << "Dummy Loop Encountered"<< std::endl;
			}
				timing(wcTimeEnd,cpuTimeEnd);	
				runtime=wcTimeEnd-wcTimeStart;				
		}
		repeat/=2;
		noFlops = N*2*repeat;
		mFlops = noFlops/runtime/1000000;
		fileoutb << N <<"\t" << mFlops << std::endl;	

		delete [] a;
		delete [] b;
		delete [] c;
		delete [] d;	
	}
	fileouta.close();
	fileoutb.close();
	std::cout << "Task 2.1 and Task 2.2 are complete" << std::endl;
	std::cout << "Running Task 2.3" << std::endl;
	// Task 2.3: Matrix Vector Multiplication
	for (int k =10; k<23; ++k)
	{
		int N = int(pow(1.5,k));
		double  **A = new double*[N];
		for (int i = 0; i<N; ++i)
			A[i] = new double[N];
	
		for (int i = 0; i<N; ++i)
		{
			for (int j = 0; j<N ; ++j)						
				A[i][j] = i*1.0+j*1.0;
		}
		double *a = new double[N];
		double *b = new double[N];
		
		for (int i = 0; i<N; ++i)
		{
			a[i] = i*1.0;
			b[i] = 0.0;
		}		
		int repeat =1;
		double runtime =0.;
		for (; runtime <0.1; repeat*=2)
		{	
			timing(wcTimeStart,cpuTimeStart);
			for (int r=0; r<repeat; ++r)
			{
				for (int i = 0; i<N; ++i)
				{
					for (int j = 0; j<N ; ++j)
						b[i] += A[i][j]*a[j];
				}
				if (b[N-1] <0.) std::cout << "Dummy Loop Encountered"<< std::endl;
			}
				timing(wcTimeEnd,cpuTimeEnd);	
				runtime=wcTimeEnd-wcTimeStart;				
		}
		repeat/=2;
		long noFlops = (2*N-1)*N*repeat;
		double mFlops = noFlops/runtime/1000000;
		
		fileoutc << N <<"\t" << mFlops << std::endl;
		delete [] a;
		delete [] b;
		for (int i = 0; i< N; ++i)
		{
			delete [] A[i];
		}
		delete [] A;
	}

	std::cout << "Task 2.3 is complete. Run the GNU Plot with plot'a-c'.dat to get the graphs " << std::endl;
	fileoutc.close();
	return 0;
}
